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Summary: We investigate the relationship between the average 
fitness decay due to single mutations and the strength of epistatic 
interactions in genetic sequences. We observe that epistatic in- 
teractions between mutations are correlated to the average fitness 
decay, both in RNA secondary structure prediction as well as in 
digital organisms replicating in silico. This correlation implies that 
during adaptation, epistasis and average mutational effect cannot 
be optimized independently. In experiments with RNA sequences 
evolving on a neutral network, the selective pressure to decrease 
the mutational load then leads to a reduction of the amount of se- 
quences with strong antagonistic interactions between deleterious 
mutations in the population. 

Keywords: epistasis, neutrality, RNA secondary structure folding, digital 
organisms 

1 Introduction 

A thorough understanding of epistatic interactions of mutations in genomes 
is becoming more and more crucial to many areas in population genetics and 
evolutionary biology. Epistasis affects linkage disequilibria (Charlesworth 
1976, Barton 1995), robustness to mutations (Lenski et al. 1999) or canal- 
ization (Nowak et al. 1997, Wagner et al. 1997, Rice 1998, Ancel & Fontana 
2000, for reviews, see Scharloo 1991, Gibson & Wagner 2000), as well as 
theories on the maintenance of sex (Kondrashov 1982, Kondrashov 1988, 
West et al. 1999). The sign of epistatic effects, that is, whether deleterious 
mutations are reinforcing (synergistic epistasis) or mitigating (antagonistic 
epistasis) also influences whether or not deleterious mutations can accumu- 
late in the genome via Muller's ratchet (Muller 1964, Felsenstein 1974). The 
consensus seems to be that synergistic epistasis can prevent the accumula- 
tion of mutations (Crow & Kimura 1978, Kondrashov 1994, but see Butcher 
1995 for a dissenting view). On the other hand, the observation of pervasive 
compensatory mutations (Moore et al. 2000), which also render the ratchet 
powerless, indicates epistasis, but not its sign. 

While the genomes of a number of organisms have been examined for 
signs of epistasis (de Visser et al. 1997a, 1997b, de Visser & Hoekstra 1998, 
Elena & Lenski 1997, Elena 1999), no general trend can be discerned except 
to say that interactions between mutations are frequent and of both signs, 
and that weak synergistic epistasis seems to prevail in eukaryotic genomes 
while viral and prokaryotic genomes show no net tendency in either direc- 
tion. Experiments to measure epistatic interactions are difficult and usually 
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yield results of weak statistical significance (West et al. 1 998). Consequently, 
even epistasis of considerable strength can conceivably be missed in vitro or 
in vivo. Here, we investigate deleterious mutations in silico, and study a fact 
which has not received much attention in the population genetics literature. 
Namely, epistasis is closely related to the geometry of the phenotype space 
(Rice 1998), which leads to interesting relations between epistasis and the 
effects of single mutations. Wagner et al. (1998) showed in a two- locus, 
two-allele model that the average effect of a single mutation is correlated to 
the degree of interaction between loci, and supported their theoretical model 
with QTL data on body weight in mice. We demonstrate here a similar rela- 
tion for a measure of epistasis which is more appropriate for high dimensional 
sequence spaces. We give both theoretical and experimental evidence that 
the strength of directional epistasis is correlated with the average deleterious 
effect of a single mutation. As a corollary of this observation, we argue that 
in situations in which there is a selective pressure to reduce the average dele- 
terious effect, this correlation leads to a reduction of the number of genomes 
with strong antagonistic effects in a population. 

2 Neutrality and Epistasis 

It is a common observation that the average fitness at a mutational distance 
n from a given reference sequence decays roughly exponentially with n (see, 
e.g., de Visser et al. 1997b, Elena & Lenski 1997, Lenski et al. 1999, and also 
our results on RNA sequences below). The simplest explanation for a perfect 
exponential is a multiplicative landscape in which every mutation diminishes 
the fitness independently by the same factor (1 — s). Here, however, we 
focus on fitness landscapes with a considerable amount of lethal mutations, 
in which case a branching process in the high-dimensional sequence space 
is a more appropriate model: if every viable sequence has a probability of 
1 — s to remain viable after suffering a single point mutation, then the mean 
fitness will decay as (neglecting fitness differences in the viable sequences) 

w{n) = (1 - s) n = e- an , (1) 

where we have defined a = — ln(l — s), and assuming that the fitness of 
our reference sequence at n = is w(0) = 1. Both explanations for the 
exponential decay have in common that the effects of subsequent mutations 
are independent, i.e., there exists no epistasis. If there are some interactions 
between mutations, we will still observe an exponential decay if antagonistic 
and synergistic interactions occur in the same proportion. If however there 
exists a bias towards either antagonistic or synergistic epistatic interactions 
(directional epistasis), then this bias will naturally appear as deviations from 



2 



the exponential decay in Eq. ([!]). While such deviations have previously been 
indicated by adding a term quadratic in n to the exponent of Eq. ([[]) (Crow 
1970, Charlesworth 1990, de Visser et al. 1997b, Elena & Lenski 1997), such 
a parameterization becomes troublesome at larger n, because w(n) could 
increase beyond the fitness of the reference sequence (for a positive coefficient 
in the quadratic term). This is avoided by the ansatz (Lenski et al. 1999) 

w (n) = e~ anP , (2) 

where (3 = 1 means there is no bias towards either form of epistatic interac- 
tions. A (3 > 1 indicates synergistic mutations prevail (mutations that are 
on average "worth" more than one independent hit), while (3 < 1 reflects 
a bias towards antagonistic mutations (mutations whose "damage" is less 
then one independent mutation) .0 Since expression @ depends only on two 
parameters, deviations from that form may arise when n grows large. 

Naively, one might assume that the decay parameter a and the epistasis 
parameter (3 are independent. Instead, we shall see that environments with 
strong selection force a trade-off between a and (3, so that one can only be 
optimized at the expense of the other. The reasoning is as follows. In a 
strongly selective environment mutations can be classified as either neutral 
or lethal, and w(n) can be thought of as the fraction of neutral sequences in 
genetic space at mutational distance n. In particular, the neutrality v of a 
sequence (the number of sequences at Hamming distance 1 with fitness 1) is 
related to the decay parameter by v = £(D — l)e~ Q , where £ is the length of 
the sequence, and D is the number of monomers. If all sequences in genetic 
space have the same u, it follows that f3 — 1. A deviation from (3=1 implies 
that some sequences have more or fewer neutral neighbors than others, giving 
rise to a correlation between a and (3. For a viable sequence with lower than 
average neutrality (higher than average a), there are comparatively fewer 
sequences close-by than there are far away, such that this sequence will have 
a small (3. Conversely, a sequence with a high neutrality (small a) will 
have comparatively more sequences close by, and (3 will be larger. We can 
make this argument more formal with a simple "conservation law", which 
only reflects that the total number of neutral sequences in genetic space is 
constant. Since for polymers of fixed length £ made from D monomers there 
are ( )(D — l) n possible n-mutants, we must have 

Y^w{n)(%D-lY = N v , (3) 

71=1 ^ ' 

1 Note that in the earlier works where a quadratic term was used, the distinction between 
the different types of directional epistasis depended on whether f3 was larger or smaller 
than 0, rather than 1. 
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where N v is the total number of neutral mutants of this wild type. Inserting 
w(n) from Eq. (|2|) yields an implicit relation between a and (3. Indeed, 
for two decay functions of two different reference sequences, with different 
parameters a, the only way in which the sum can yield the quantity N v 
(which is independent of the respective reference sequence) is that the two 
decay functions must have different parameters (3 as well. However, the 
implicit relation depends on the ansatz Eq. (§) being correct for all n, which 
is not necessarily the case. Alternatively, we may consider only sequences 
with up to d mutations, in which case we can write 



where N v% d is the number of neutral sequences in a sphere of radius d around 
the reference sequence with w(0) = 1. N V) d depends on the particular refer- 
ence sequence chosen, and therefore cannot be regarded a constant. However, 
for d not too small, we may replace it by its average (N u ^) over all viable 
reference sequences. If we take into account a sufficiently large region of geno- 
type space, we should find roughly the same number of neutral mutants for 
each viable sequence, since, as d approaches £, the quantity N„ t d approaches 
N u , which in turn is independent of the reference sequence. Hence, Eq. ([|) 
predicts a similar relation between a and f3, and the two different predictions 
approach each other as d — > I. Predictions based on Eqs. (|3]) and @ are used 
below to compare to our empirical results. 

Although the above argument strictly holds only under the assumption 
that mutations are either neutral or lethal, it is not unreasonable to assume 
that a similar (possibly weaker) correlation between a and (3 exists also in 
more general cases, where slightly deleterious or even advantageous muta- 
tions are possible. In that case, under the presence of epistasis, there will 
still be regions in genotype space in which the number of less-deleterious 
mutations is higher, and other regions in which it is lower than average. The 
decay function w(n) of a sequence from a region that is rich in non-lethal mu- 
tations would have a lower a, but would be inevitably more synergistic than 
the decay function of a sequence from a region poor in non-lethal mutations. 
Our results with digital organisms (see below) support this reasoning. 

3 Experimental Evidence 

Accurate data for the decay parameter a and the epistasis parameter j3 for 
biological organisms are rare, which makes our hypothesis difficult to test. 
A few well-studied systems have emerged which are accessible in silico. 
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We studied RNA secondary structure prediction using the Vienna RNA 
package, version 1.3.1., with the default setup (Hofacker et al. 1994). We 
calculated the decay of the average number of neutral folds as a function of 
the Hamming distance for 100 random RNA sequences of length £ = 76. The 
parameters a and (3 were determined as follows. We obtained a exactly from 
the fraction of neutral one-mutants. In addition, we sampled the function 
w(n) for Hamming distances up to n = 8, by calculating the structure of up 
to 10 6 random neighbors of the required Hamming distance. The quantity 
(3 was then determined from a nonlinear fit of —arfi to the logarithm of 
w(n). A plot of (3 versus a (Fig. [I]) shows a significant correlation, with a 
correlation coefficient of r = —0.817 (p < 0.01). 

According to Eqs. (0) and (|]), we can predict the relationship between 
a and (3 if we compare the decay functions of sequences that are mutually 
neutral. For RNA folding, this means we have to determine a and (3 for a set 
of sequences that fold into the same structure. We performed experiments 
with the RNA sequences of length I = 18 used by van Nimwegen et al. 
(1999) (for these sequences we altered the default setup by setting the free 
energies of dangling ends to zero.) For the particular case that all bonds 
are of the purine- pyrimidine type (G-C, G-U, A-U), two separate neutral 
networks (a neutral network is a network of neutral genotypes connected 
to each other by one-point mutations) were found by van Nimwegen et al, 
consisting of 51,028 and 5,169 sequences, respectively. For each such set of 
neutral sequences, Eq. (|3|) predicts the correlation without free parameter, 
as long as the number of all neutral sequences N v is known. In order to 
estimate N u , we generated 10 8 random sequences of length t = 18, of which 
10,961 sequences folded correctly. From this, we estimated N v = 7.5 x 10 6 
neutral sequences out of the total 6.9 x 10 10 sequences of length £ = 18. 
Using this number, Eq. (|3|) predicts the solid line in Fig. which describes 
the correlation well. The approach based on averaging over local neutral 
sequences around the reference sequence [Eq. (|])] gives rise to the dotted 
line in Fig. |2| and shows that this too predicts the correlation fairly well (for 
this second approach, we used d = 8, and from Monte Carlo simulations for 
1000 reference sequences, we obtained (N Ut8 ) = 2.4 x 10 5 ). 

As our second test case, we analyzed the correlation between a and (3 in 
digital organisms (Adami 1998, Adami et al. 2000). Digital organisms are 
self-replicating computer programs that mutate and evolve. Their fitness is 
determined by the ratio of their CPU speed and their gestation time. The lat- 
ter is given by the number of instructions that have to be executed in order to 
generate a fully functional offspring. The CPU speed increases when the dig- 
ital organisms perform logical operations on numbers that they can acquire 
from their environment. The gestation time, on the other hand, increases 
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either if the organisms employ less efficient mechansims of self-replication, 
or if they accumulate instructions that are involved in the completion of the 
above mentioned logical operations. The digital organisms with the highest 
fitness thus have a very efficient copy mechanism, and perform a large num- 
ber of logical operations with a comparatively small number of additional 
instructions. Lenski et al. (1999) measured the decay of the mean fitness as 
a function of the number of mutations accumulated in such digital organisms, 
and obtained a and (3 from a fit of Eq. (|2]) to the measured decay functions. 
They studied 174 different genomes consisting of two groups of 87 genomes 
each. The first group of organisms evolved in 87 independent experiments 
in a complex environment, while the second group was obtained by transfer- 
ing these organisms to an environment which favoured simple genomes, and 
allowing the organisms to adapt to this more simple environment. 

A statistical analysis revealed a significant correlation between the decay 
parameter a and the parameter of directional epistasis f3 for both the complex 
and the simple digital organisms (Table |l|). However, in addition to said 
correlation, we found a correlation between the a and both the genome length 
£ and the log fitness for complex and simple organisms, as well as a correlation 
between (3 and £ in the case of the complex organisms. Hence, for the complex 
organisms, we cannot rule out the possibility that the correlation between 
a and (3 merely reflects an underlying correlation of both quantities with 
length. In the case of the simple organisms, where we do not see a correlation 
between (3 and £, we can assume that the correlation between a and (3 is 
genuine. To provide further evidence, we examined a reduced data set of all 
48 simple organisms with a length between £ = 14 and £ = 16. These 48 
organisms were also of comparable fitness. In that data set, we found an even 
stronger correlation between a and (3, while the correlation between either 
of the two quantities and length or fitness was insignificant (Table [I]). It was 
not possible to study a similar reduced data set for the complex organisms, 
because the variations in length were too large (the length varied between 
£ = 20 and £ = 314 among the 87 genomes). Although the data from the 
digital organisms is potentially biased because our reference sequences were 
evolved rather than random as in the RNA case, we believe this bias to be 
insubstantial in at least the reduced data set, which contains only organisms 
with very similar phenotypes. Hence, just as for our RNA data, the data 
from the digital organisms supports our hypothesis of a genuine correlation 
between a and (3. 
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4 Adaptation of Epistasis through Correlated 
Response 

The correlation between neutrality and epistasis implies that if one of them 
is subject to selective pressures, the other will be affected as well due to cor- 
related response. Van Nimwegen et al. (f999) have shown that a population 
evolving on a neutral network reduces its genetic load by moving into the 
regions of high neutrality in sequence space. In particular, given a random 
population of molecules on such a network, evolution tends to increase the 
neutrality in the population, effectively pushing the population into the cen- 
ter of the neutral network. Because of the correlation between neutrality 
and epistasis, we expect this dynamic to lead to a reduction of antagonis- 
tic epistatic effects. To verify this hypothesis, we carried out evolutionary 
experiments with the RNA sequences of length 18 from the previous section. 

We performed one flow-reactor run for each of the two networks found by 
van Nimwegen et al, starting with an initial population of 1,000 sequences 
chosen at random from the respective network. We set the replication rates 
such that sequences folding into the target structure replicated on average 
once per unit time, while the replication rate of all other sequences was set to 
10~ 6 per unit time. All sequences had a probability of \x = 0.5 to suffer one 
random point mutation per replication event. The possibility of several point 
mutations per replication event was eliminated, to guarantee that the popu- 
lation could not leave the specified neutral network. The epistasis parameter 
(3 was determined for every sequence in the population every two hundred 
generations, while the population neutrality was monitored constantly. The 
population neutrality v is the average neutrality of all sequences currently 
in the population. In Fig. |3|, we present the results from the run on the 
larger of the two networks. The neutrality of the initial population coincides 
with the network neutrality (the average neutrality of all sequences on the 
network), which is to be expected for a random initial population. Over the 
course of evolution, the average population neutrality rose to the predicted 
equilibrium value (given by the spectral radius of the connectivity matrix, 
see van Nimwegen et al. 1999.) As expected, the average epistasis param- 
eter (3 increased significantly as well. Results on the second network were 
qualitatively identical, with f3 increasing from 0.78 to around 0.86. Thus, 
antagonistic epistasis is reduced during adaptation for reduced mutational 
load on a neutral network. 

The above reasoning depends of course on the assumption that a popula- 
tion remains on a single neutral network. In a more realistic scenario, where 
peaks of different heights are present, the main effect of selection will be to 
increase the fitness, rather than to increase a. However, once a local opti- 
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mum has been reached, we can expect the dynamic described above to take 
place, as long as there exists some neutrality at the local peak. If there was 
a correlation between the fitness and a or f3, then we would see additionally 
a correlated response in a or (3 as the fitness is being maximized. Neverthe- 
less, while the correlation between a and (3 reported here is a general result 
that follows from geometric constraints on the landscape, no such constraints 
exist between the fitness and a or (3. The parameters a and (3 measure the 
amount and the distribution of neutral or nearly neutral sequences in the 
neighborhood of a reference sequence. There is no reason why sequences 
with high fitness should generally be found in regions with particularly high 
or low neutrality, or with a particular distribution of the neutral sequences. 
This is not to say, however, that such a correlation cannot exist in special 
cases (see, e.g., our data from the complex digital organisms, Table [Q where 
a is correlated with fitness, but f3 is not). 

5 Conclusions 

Epistasis plays an important role in evolutionary theory, but remains em- 
pirically largely unexplored. Using secondary structure prediction of RNA 
sequences as well as digital organisms evolving in silico, we have demon- 
strated a correlation between two important parameters of realistic genetic 
fitness landscapes: the average deleterious effect of single mutations and the 
strength of directional epistasis. In conjunction with the results from Wagner 
et al. (1998) for mice, and the two different theoretical explanations (two- 
locus model in the case of Wagner et al, sequence space based model here), 
we can expect this correlation to be an ubiquitious phenomenon, present in 
many natural and artificial fitness landscapes. 

The correlation, coupled with the selective pressure which forces random 
sequences in a neutral network to cluster in the dense areas of the network, 
leads to a reduction of strong antagonistic epistasis in a population. The 
nature of this result is purely geometric: as a population tries to reduce 
the average effect of single mutations, the effect of multiple mutations is 
inevitably worsened as long as there exist some inhomogeneities in the effect 
of single mutations across the genotype space. As the result of this geometric 
constraint, a member of an evolved population will have, on average, a higher 
(3 than a random sample of the fitness landscape would indicate. 

It is well known that antagonistic epistasis favors the accumulation of 
deleterious mutations as well as the operation of Muller's ratchet. Since in 
such a situation sexual recombination (within a fixed environment) tends 
to worsen the loss of information, recombination is unlikely to evolve or be 
maintained. The mechanism described here may thus provide a path towards 
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an environment more conducive to the evolution of recombination. 
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Table 1: The correlation r and p- value between decay parameter a, epistasis 
parameter (3, length £ and the logarithm of the fitness In w in the data from 
(Lenski et al. 1999). The "Complex" and the "Simple" data set consist each 
of 87 digital organisms, the "Reduced" data set consists of all 48 organisms 
of length between 14 and 16 taken from the "Simple" data set. 
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Figure 1: Epistasis parameter (3 versus the decay parameter a in random 
RNA sequences of length £ = 76. 
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Figure 2: Epistasis parameter (3 versus the decay parameter a for sequences 
on the larger of the two neutral networks from (van Nimwegen et al. 1999). 
The solid line represents the prediction that follows from Eq. (|3[), and the 
dotted line represents the prediction from Eq. (f|), with d — 8. 
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Figure 3: Evolution of neutrality and epistasis as a function of time (in gen- 
erations). The lower graph shows the convergence of the population neutral- 
ity to the value predicted by the spectral radius of the connectivity matrix. 
The upper graph shows the change of (3, averaged over the population, in the 
same run. For all data points in the upper graph, the standard error of the 
mean does not exceed the size of the symbols. 
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